{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
       "@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
       "@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
       "@import url('http://fonts.googleapis.com/css?family=Fira_sans');\n",
       "\n",
       "    div.cell{\n",
       "        width: 900px;\n",
       "        margin-left: 0% !important;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    div.text_cell code {\n",
       "        background: transparent;\n",
       "        color: #000000;\n",
       "        font-weight: 600;\n",
       "        font-size: 11pt;\n",
       "        font-style: bold;\n",
       "        font-family:  'Source Code Pro', Consolas, monocco, monospace;\n",
       "   }\n",
       "    h1 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "\t}\n",
       "\t\n",
       "    div.input_area {\n",
       "        background: #F6F6F9;\n",
       "        border: 1px solid #586e75;\n",
       "    }\n",
       "\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 30pt;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h2 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "        text-align: left;\n",
       "    }\n",
       "    .text_cell_render h2 {\n",
       "        font-weight: 200;\n",
       "        font-size: 16pt;\n",
       "        font-style: italic;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1.5em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h3 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h3 {\n",
       "        font-weight: 200;\n",
       "        font-size: 14pt;\n",
       "        line-height: 100%;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 2em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    }\n",
       "    h4 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h4 {\n",
       "        font-weight: 100;\n",
       "        font-size: 14pt;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    h5 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 200;\n",
       "        font-style: normal;\n",
       "        color: #1d3b84;\n",
       "        font-size: 16pt;\n",
       "        margin-bottom: 0em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: 'Fira sans', verdana,arial,sans-serif;\n",
       "        line-height: 125%;\n",
       "        font-size: 115%;\n",
       "        text-align:justify;\n",
       "        text-justify:inter-word;\n",
       "    }\n",
       "    div.output_subarea.output_text.output_pyout {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_subarea.output_stream.output_stdout.output_text {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_wrapper{\n",
       "        margin-top:0.2em;\n",
       "        margin-bottom:0.2em;\n",
       "}\n",
       "\n",
       "    code{\n",
       "      font-size: 70%;\n",
       "    }\n",
       "    .rendered_html code{\n",
       "    background-color: transparent;\n",
       "    }\n",
       "    ul{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li li{\n",
       "        padding-left: 0.2em; \n",
       "        margin-bottom: 0.2em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    ol{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ol li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    a:link{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:visited{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:hover{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:focus{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:active{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    .rendered_html :link {\n",
       "       text-decoration: underline; \n",
       "    }\n",
       "    .rendered_html :hover {\n",
       "       text-decoration: none; \n",
       "    }\n",
       "    .rendered_html :visited {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :focus {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :active {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "    } \n",
       "    hr {\n",
       "      color: #f3f3f3;\n",
       "      background-color: #f3f3f3;\n",
       "      height: 1px;\n",
       "    }\n",
       "    blockquote{\n",
       "      display:block;\n",
       "      background: #fcfcfc;\n",
       "      border-left: 5px solid #c76c0c;\n",
       "      font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "      width:680px;\n",
       "      padding: 10px 10px 10px 10px;\n",
       "      text-align:justify;\n",
       "      text-justify:inter-word;\n",
       "      }\n",
       "      blockquote p {\n",
       "        margin-bottom: 0;\n",
       "        line-height: 125%;\n",
       "        font-size: 100%;\n",
       "      }\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    scale:100,\n",
       "                        availableFonts: [],\n",
       "                        preferredFont:null,\n",
       "                        webFont: \"TeX\",\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# format the book\n",
    "%matplotlib inline\n",
    "import sys\n",
    "from __future__ import division, print_function\n",
    "import sys\n",
    "sys.path.insert(0,'../code')\n",
    "import book_format\n",
    "book_format.load_style('../code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Hierarchical Model\n",
    "\n",
    "## The Geiger counter problem\n",
    "\n",
    "I got the idea for the following problem from Tom Campbell-Ricketts,\n",
    "author of the Maximum Entropy blog at\n",
    "<http://maximum-entropy-blog.blogspot.com>. And he got the idea from\n",
    "E. T. Jaynes, author of the classic `Probability Theory: The Logic\n",
    "of Science`:\n",
    "\n",
    "> Suppose that a radioactive source emits particles toward a Geiger\n",
    "> counter at an average rate of $r$ particles per second, but the\n",
    "> counter only registers a fraction, $f$, of the particles that hit it.\n",
    "> If $f$ is 10% and the counter registers 15 particles in a one second\n",
    "> interval, what is the posterior distribution of $n$, the actual number\n",
    "> of particles that hit the counter, and $r$, the average rate particles\n",
    "> are emitted?\n",
    "\n",
    "To get started on a problem like this, think about the chain of\n",
    "causation that starts with the parameters of the system and ends with\n",
    "the observed data:\n",
    "\n",
    "1.  The source emits particles at an average rate, $r$.\n",
    "\n",
    "2.  During any given second, the source emits $n$ particles toward the\n",
    "    counter.\n",
    "\n",
    "3.  Out of those $n$ particles, some number, $k$, get counted.\n",
    "\n",
    "The probability that an atom decays is the same at any point in time, so\n",
    "radioactive decay is well modeled by a Poisson process. Given $r$, the\n",
    "distribution of $n$ is Poisson distribution with parameter $r$.\n",
    "\n",
    "And if we assume that the probability of detection for each particle is\n",
    "independent of the others, the distribution of $k$ is the binomial\n",
    "distribution with parameters $n$ and $f$.\n",
    "\n",
    "Given the parameters of the system, we can find the distribution of the\n",
    "data. So we can solve what is called the **forward problem**.\n",
    "\n",
    "Now we want to go the other way: given the data, we want the\n",
    "distribution of the parameters. This is called the **inverse problem**. \n",
    "And if you can solve the forward problem, you can use\n",
    "Bayesian methods to solve the inverse problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Start simple\n",
    "\n",
    "![Posterior distribution of $n$ for three values of\n",
    "$r$.](figs/jaynes1.pdf)\n",
    "\n",
    "[fig.jaynes1]\n",
    "\n",
    "Let’s start with a simple version of the problem where we know the value\n",
    "of $r$. We are given the value of $f$, so all we have to do is estimate\n",
    "$n$.\n",
    "\n",
    "I define a Suite called `Detector` that models the behavior\n",
    "of the detector and estimates $n$.\n",
    "\n",
    "```python\n",
    "class Detector(thinkbayes.Suite):\n",
    "\n",
    "    def __init__(self, r, f, high=500, step=1):\n",
    "        pmf = thinkbayes.MakePoissonPmf(r, high, step=step)\n",
    "        thinkbayes.Suite.__init__(self, pmf, name=r)\n",
    "        self.r = r\n",
    "        self.f = f\n",
    "```\n",
    "\n",
    "If the average emission rate is $r$ particles per second, the\n",
    "distribution of $n$ is Poisson with parameter $r$. `high` and\n",
    "`step` determine the upper bound for $n$ and the step size\n",
    "between hypothetical values.\n",
    "\n",
    "Now we need a likelihood function:\n",
    "\n",
    "```python\n",
    "# class Detector\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        k = data\n",
    "        n = hypo\n",
    "        p = self.f\n",
    "\n",
    "        return thinkbayes.EvalBinomialPmf(k, n, p)\n",
    "```\n",
    "`data` is the number of particles detected, and\n",
    "`hypo` is the hypothetical number of particles emitted, $n$.\n",
    "\n",
    "If there are actually $n$ particles, and the probability of detecting\n",
    "any one of them is $f$, the probability of detecting $k$ particles is\n",
    "given by the binomial distribution.\n",
    "\n",
    "That’s it for the Detector. We can try it out for a range of values of\n",
    "$r$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "OverflowError",
     "evalue": "int too large to convert to float",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mOverflowError\u001b[0m                             Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-d79750c08982>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mr\u001b[0m \u001b[1;32min\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m100\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m250\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m400\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 8\u001b[1;33m     \u001b[0msuite\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mDetector\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstep\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      9\u001b[0m     \u001b[0msuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mUpdate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     10\u001b[0m     \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMaximumLikelihood\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Users\\rlabbe\\Dropbox\\ThinkBayes\\code\\jaynes.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, r, f, high, step)\u001b[0m\n\u001b[0;32m    132\u001b[0m         \u001b[0mstep\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mstep\u001b[0m \u001b[0msize\u001b[0m \u001b[0mbetween\u001b[0m \u001b[0mhypothetical\u001b[0m \u001b[0mvalues\u001b[0m \u001b[0mof\u001b[0m \u001b[0mn\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    133\u001b[0m         \"\"\"\n\u001b[1;32m--> 134\u001b[1;33m         \u001b[0mpmf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mthinkbayes\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMakePoissonPmf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhigh\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstep\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mstep\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    135\u001b[0m         \u001b[0mthinkbayes\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__init__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpmf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    136\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mr\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Users\\rlabbe\\Dropbox\\ThinkBayes\\code\\thinkbayes.py\u001b[0m in \u001b[0;36mMakePoissonPmf\u001b[1;34m(lam, high, step)\u001b[0m\n\u001b[0;32m   1563\u001b[0m     \u001b[0mpmf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mPmf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1564\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhigh\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstep\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1565\u001b[1;33m         \u001b[0mp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mEvalPoissonPmf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlam\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1566\u001b[0m         \u001b[0mpmf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSet\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1567\u001b[0m     \u001b[0mpmf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mNormalize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Users\\rlabbe\\Dropbox\\ThinkBayes\\code\\thinkbayes.py\u001b[0m in \u001b[0;36mEvalPoissonPmf\u001b[1;34m(k, lam)\u001b[0m\n\u001b[0;32m   1550\u001b[0m     \u001b[1;31m# return scipy.stats.poisson.pmf(k, lam)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1551\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1552\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mlam\u001b[0m \u001b[1;33m**\u001b[0m \u001b[0mk\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mmath\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mlam\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m/\u001b[0m \u001b[0mmath\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfactorial\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1553\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1554\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mOverflowError\u001b[0m: int too large to convert to float"
     ]
    }
   ],
   "source": [
    "from jaynes import Detector\n",
    "\n",
    "f = 0.1\n",
    "k = 15\n",
    "\n",
    "\n",
    "for r in [100, 250, 400]:\n",
    "    suite = Detector(r, f, step=1)\n",
    "    suite.Update(k)\n",
    "    print(suite.MaximumLikelihood())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure [fig.jaynes1] shows the posterior distribution of $n$ for several\n",
    "given values of $r$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make it hierarchical\n",
    "\n",
    "In the previous section, we assume $r$ is known. Now let’s relax that\n",
    "assumption. I define another Suite, called `Emitter`, that\n",
    "models the behavior of the emitter and estimates $r$:\n",
    "\n",
    "```python\n",
    "class Emitter(thinkbayes.Suite):\n",
    "\n",
    "    def __init__(self, rs, f=0.1):\n",
    "        detectors = [Detector(r, f) for r in rs]\n",
    "        thinkbayes.Suite.__init__(self, detectors)\n",
    "```\n",
    "\n",
    "`rs` is a sequence of hypothetical value for $r$.\n",
    "`detectors` is a sequence of Detector objects, one for each\n",
    "value of $r$. The values in the Suite are Detectors, so Emitter is a\n",
    "`**meta-Suite**`; that is, a Suite that contains other Suites\n",
    "as values.\n",
    "\n",
    "To update the Emitter, we have to compute the likelihood of the data\n",
    "under each hypothetical value of $r$. But each value of $r$ is\n",
    "represented by a Detector that contains a range of values for $n$.\n",
    "\n",
    "To compute the likelihood of the data for a given Detector, we loop\n",
    "through the values of $n$ and add up the total probability of $k$.\n",
    "That’s what `SuiteLikelihood` does:\n",
    "\n",
    "```python\n",
    "# class Detector\n",
    "\n",
    "    def SuiteLikelihood(self, data):\n",
    "        total = 0\n",
    "        for hypo, prob in self.Items():\n",
    "            like = self.Likelihood(data, hypo)\n",
    "            total += prob * like\n",
    "        return total\n",
    "```\n",
    "\n",
    "Now we can write the Likelihood function for the Emitter:\n",
    "\n",
    "```python\n",
    "# class Detector\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        detector = hypo\n",
    "        like = detector.SuiteLikelihood(data)\n",
    "        return like\n",
    "```\n",
    "\n",
    "Each `hypo` is a Detector, so we can invoke\n",
    "`SuiteLikelihood` to get the likelihood of the data under the\n",
    "hypothesis.\n",
    "\n",
    "After we update the Emitter, we have to update each of the Detectors,\n",
    "too.\n",
    "\n",
    "```python\n",
    "# class Detector\n",
    "\n",
    "    def Update(self, data):\n",
    "        thinkbayes.Suite.Update(self, data)\n",
    "\n",
    "        for detector in self.Values():\n",
    "            detector.Update()\n",
    "```\n",
    "\n",
    "A model like this, with multiple levels of Suites, is called\n",
    "**hierarchical**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A little optimization\n",
    "\n",
    "You might recognize `SuiteLikelihood`; we saw it in\n",
    "Section [suitelike]. At the time, I pointed out that we didn’t really\n",
    "need it, because the total probability computed by\n",
    "`SuiteLikelihood` is exactly the normalizing constant\n",
    "computed and returned by `Update`.\n",
    "\n",
    "So instead of updating the Emitter and then updating the Detectors, we\n",
    "can do both steps at the same time, using the result from\n",
    "`Detector.Update` as the likelihood of Emitter.\n",
    "\n",
    "Here’s the streamlined version of `Emitter.Likelihood`:\n",
    "\n",
    "```python\n",
    "# class Emitter\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        return hypo.Update(data)\n",
    "```\n",
    "\n",
    "And with this version of `Likelihood` we can use the default\n",
    "version of `Update`. So this version has fewer lines of code,\n",
    "and it runs faster because it does not compute the normalizing constant\n",
    "twice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting the posteriors\n",
    "\n",
    "![Posterior distributions of $n$ and $r$.](figs/jaynes2.pdf)\n",
    "\n",
    "[fig.jaynes2]\n",
    "\n",
    "After we update the Emitter, we can get the posterior distribution of\n",
    "$r$ by looping through the Detectors and their probabilities:\n",
    "\n",
    "```python\n",
    "# class Emitter\n",
    "\n",
    "    def DistOfR(self):\n",
    "        items = [(detector.r, prob) for detector, prob in self.Items()]\n",
    "        return thinkbayes.MakePmfFromItems(items)\n",
    "```\n",
    "\n",
    "`items` is a list of values of $r$ and their probabilities.\n",
    "The result is the Pmf of $r$.\n",
    "\n",
    "To get the posterior distribution of $n$, we have to compute the mixture\n",
    "of the Detectors. We can use `thinkbayes.MakeMixture`, which\n",
    "takes a meta-Pmf that maps from each distribution to its probability.\n",
    "And that’s exactly what the Emitter is:\n",
    "\n",
    "```python\n",
    "# class Emitter\n",
    "\n",
    "    def DistOfN(self):\n",
    "        return thinkbayes.MakeMixture(self)\n",
    "```\n",
    "\n",
    "Figure [fig.jaynes2] shows the results. Not surprisingly, the most\n",
    "likely value for $n$ is 150. Given $f$ and $n$, the expected count is\n",
    "$k = f n$, so given $f$ and $k$, the expected value of $n$ is $k / f$,\n",
    "which is 150.\n",
    "\n",
    "And if 150 particles are emitted in one second, the most likely value of\n",
    "$r$ is 150 particles per second. So the posterior distribution of $r$ is\n",
    "also centered on 150.\n",
    "\n",
    "The posterior distributions of $r$ and $n$ are similar; the only\n",
    "difference is that we are slightly less certain about $n$. In general,\n",
    "we can be more certain about the long-range emission rate, $r$, than\n",
    "about the number of particles emitted in any particular second, $n$.\n",
    "\n",
    "You can download the code in this chapter from\n",
    "<http://thinkbayes.com/jaynes.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "The Geiger counter problem demonstrates the connection between causation\n",
    "and hierarchical modeling. In the example, the emission rate $r$ has a\n",
    "causal effect on the number of particles, $n$, which has a causal effect\n",
    "on the particle count, $k$.\n",
    "\n",
    "The hierarchical model reflects the structure of the system, with causes\n",
    "at the top and effects at the bottom.\n",
    "\n",
    "1.  At the top level, we start with a range of hypothetical values for\n",
    "    $r$.\n",
    "\n",
    "2.  For each value of $r$, we have a range of values for $n$, and the\n",
    "    prior distribution of $n$ depends on $r$.\n",
    "\n",
    "3.  When we update the model, we go bottom-up. We compute a posterior\n",
    "    distribution of $n$ for each value of $r$, then compute the\n",
    "    posterior distribution of $r$.\n",
    "\n",
    "So causal information flows down the hierarchy, and inference flows up."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "This exercise is also inspired by an example in Jaynes,\n",
    "*Probability Theory*.\n",
    "\n",
    "Suppose you buy a mosquito trap that is supposed to reduce the\n",
    "population of mosquitoes near your house. Each week, you empty the trap\n",
    "and count the number of mosquitoes captured. After the first week, you\n",
    "count 30 mosquitoes. After the second week, you count 20 mosquitoes.\n",
    "Estimate the percentage change in the number of mosquitoes in your yard.\n",
    "\n",
    "To answer this question, you have to make some modeling decisions. Here\n",
    "are some suggestions:\n",
    "\n",
    "-   Suppose that each week a large number of mosquitoes, $N$, is bred in\n",
    "    a wetland near your home.\n",
    "\n",
    "-   During the week, some fraction of them, $f_1$, wander into your\n",
    "    yard, and of those some fraction, $f_2$, are caught in the trap.\n",
    "\n",
    "-   Your solution should take into account your prior belief about how\n",
    "    much $N$ is likely to change from one week to the next. You can do\n",
    "    that by adding a level to the hierarchy to model the percent change\n",
    "    in $N$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
